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Abstract 

In  this  paper,  we  present  a  comprehensive  non-isothermal,  one-dimensional  model  of  the  cathode  side  of  a  Polymer  Electrolyte  Fuel  Cell.  We 
explicitly  include  the  catalyst  layer,  gas  diffusion  layer  and  the  membrane.  The  catalyst  layer  and  gas  diffusion  layer  are  characterized  by  several 
measurable  microstructural  parameters.  We  model  all  three  phases  of  water,  with  a  view  to  capturing  the  effect  that  each  has  on  the  performance  of 
the  cell.  A  comparison  with  experiment  is  presented,  demonstrating  excellent  agreement,  particularly  with  regard  to  the  effects  of  water  activity  in 
the  channels  and  how  it  impacts  flooding  and  membrane  hydration.  We  present  several  results  pertaining  to  the  effects  of  water  on  the  current  density 
(or  cell  voltage),  demonstrating  the  role  of  micro-structure,  liquid  water  removal  from  the  channel,  water  activity,  membrane  and  gas  diffusion  layer 
thickness  and  channel  temperature.  These  results  provide  an  indication  of  the  changes  that  are  required  to  achieve  optimal  performance  through 
improved  water  management  and  MEA-component  design.  Moreover,  with  its  level  of  detail,  the  model  we  develop  forms  an  excellent  basis  for  a 
multi-dimensional  model  of  the  entire  membrane  electrode  assembly. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

The  fuel  cell,  in  its  various  forms,  is  a  promising  alternative  to 
the  use  of  hydrocarbon  fuels  for  energy  production,  motivated 
largely  by  concern  for  the  environment.  The  proton  exchange 
membrane  fuel  cell  (PEMFC)  receives  a  great  deal  of  attention, 
particularly  because  it  attains  a  relatively  high  power  density  and 
electrical  efficiency.  Before  the  benefits  of  this  technology  can  be 
realized,  several  important  practical  issues  need  to  be  resolved, 
with  the  overarching  aims  of  reducing  the  manufacturing  cost 
and  improving  the  performance  and  durability.  A  key  concern, 
which  we  address  in  this  paper,  is  the  management  of  water  in 
its  various  phases. 

The  heart  of  a  single  PEMFC  is  composed  of  a  proton  con¬ 
ducting  membrane  sandwiched  between  an  anode  and  cathode. 
The  anode  and  cathode  are  thin  sheets  of  porous  graphite  paper, 
onto  which  a  catalyst  layer,  composed  roughly  of  carbon-grain 
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supported  Platinum  particles,  is  pressed.  The  basic  principle  of 
the  PEMFC,  or  indeed  any  other  fuel  cell,  is  to  convert  chemi¬ 
cal  into  electrical  energy.  At  the  anode,  hydrogen  is  fed  through 
the  graphite  paper,  termed  the  gas  diffusion  layer  (GDL),  and 
is  oxidized  in  the  anode  catalyst  layer  (ACL),  producing  hy¬ 
drogen  ions  and  free  electrons.  The  electrons  pass  through  an 
external  circuit  and  the  ions  through  the  membrane,  eventually 
combining  with  electrons  and  oxygen  in  the  cathode  catalyst 
layer  (CCL)  to  form  water.  The  oxygen  is  fed  through  a  GDL 
at  the  cathode  side.  The  membrane  electrode  assembly  (MEA) 
is  depicted  in  Fig.  1.  Because  of  the  low  operating  temperatures 
the  Platinum  (Pt)  catalyst  is  required  to  enhance  the  rate  of  re¬ 
action  by  providing  an  alternative  pathway  to  the  final  products. 
The  cost  of  Pt  is  clearly  a  concern;  though  the  required  amount 
of  Pt  has  been  reduced  ten-fold  reduction  in  the  last  few  years, 
the  percentage  that  is  utilized  is  still  very  low  (approximately 
20%). 

The  catalyst  layers  are  amongst  the  most  complex  compo¬ 
nents  of  the  fuel  cell.  Through  the  solid  components  (carbon 
grains,  supporting  smaller  catalyst  particles,  and  electrolyte), 
flow  components  of  air  and  water  vapor.  Liquid  water  resides 
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Fig.  1 .  Schematic  representation  of  the  MEA  and  a  description  of  the  geometry. 


in  the  gas  pores  and  dissolved  water  within  the  electrolyte.  The 
electrons  migrate  through  the  carbon  components  of  the  layer 
while  the  protons  migrate  through  pathways  provided  by  the 
electrolyte.  The  latter  is  typically  Nafion  (a  sulphonated  Teflon), 
which  contains  both  hydrophobic  and  hydrophilic  regions,  the 
second  of  which  are  vital  in  ensuring  the  free  flow  of  the  protons. 
The  CCL  is  the  more  studied  of  the  two  catalyst  layers  for  several 
reasons.  The  oxygen-reduction  reaction  (ORR)  takes  place  at  a 
much  slower  rate  than  the  hydrogen  oxidation  at  the  anode  and 
produces  water.  The  hydrophilic  regions  in  the  CCL  can  cause 
water  retention,  which  restricts  the  ingress  of  oxygen.  Moreover, 
reaction  is  limited  by  the  availability  of  Pt,  and  can  only  occur 
at  points  of  contact  between  the  Pt,  carbon  and  electrolyte.  Be¬ 
cause  of  the  importance  of  the  CCL  it  is  included  explicitly  in 
our  model,  with  an  attempt  to  capture  some  of  the  micro-scale 
features  of  its  structure. 

The  properties  of  the  membrane  and  the  GDL  also  play  a 
central  role  in  the  control  of  water.  It  is  essential  that  the  protons 
generated  at  the  anode  can  reach  the  cathode  freely,  and  in  order 
to  attain  the  high  conductivity  necessary  to  achieve  this,  the 
membrane  must  remain  well  hydrated.  The  GDL  must  allow 
sufficient  oxygen  from  the  cathode  channel  to  enter  the  CCL, 
but  must  be  sufficiently  hydrophobic  to  expel  any  build-up  of 
liquid  water. 

The  preceding  discussion  underlines  the  unquestionable  im¬ 
portance  of  water  management,  which  in  turn  explains  the  grow¬ 
ing  emphasis  placed  on  (particularly)  liquid  water  in  modelling 
studies,  [1-9].  Though  many  of  the  issues  are  difficult  to  sepa¬ 
rate,  they  may  be  summarized  as  follows: 

1 .  Flooding  of  the  CCL  and  GDL  and  its  impact  on  the  perfor¬ 
mance; 

2.  The  relationship  between  humidification  of  the  cell  and  per¬ 
formance; 

3.  Hydration  of  the  membrane  and  electrolyte  (ionomer); 

4.  The  trend  of  PEMFC  operating  conditions  toward  higher  tem¬ 
perature  and  lower  water  activity,  and  their  impact  on  satu¬ 
ration  levels  (optimal  conditions  for  performance). 


One  of  the  difficulties  associated  with  PEMFC  modelling  lies 
in  accurately  describing  the  structure  of  the  CCL.  The  closest 
approximation  to  the  structure  is  one  of  clumps  of  carbon  grains 
(agglomerates)  coated  with  and  connected  by  electrolyte,  [6,10- 
14].  A  widely  employed  approach  in  modelling  the  CCL  is  to 
assume  that  the  agglomerates  are  spherical  in  shape,  [6,7, 15-20], 
and  to  incorporate  the  agglomerate-level  activity  into  a  homo¬ 
geneous  model  by  assuming  that  a  thin  film  of  electrolyte  (or 
water)  introduces  a  local  resistance  to  the  oxygen — motivated 
by  the  well-established  theory  of  porous  catalysts  (see  Chapter 
3  of  [21]).  The  three  main  variants  are  as  follows: 

•  Include  external  resistance  to  the  oxygen  movement  arising 
from  the  film(s) — this  yields  an  expression  for  the  limiting 
current  density  (or  reaction  rate)  as  a  function  of  the  agglom¬ 
erate  characteristics,  [3,6,7,15,18]. 

•  Include  internal  resistance  to  the  oxygen  movement,  due 
to  flooding  of  the  agglomerates  with  electrolyte  or  liquid 
water — this  yields  an  effectiveness  factor  for  the  oxygen  dif¬ 
fusion  coefficient,  [16,17], 

•  Include  both  internal  and  external  resistances,  as  described 
above,  [19,20], 

Explicit  incorporation  of  such  resistances  distinguishes  these 
models  from  those  in  [22-28],  which  treat  the  CCL  as  a  homo¬ 
geneous  porous  medium,  i.e. ,  each  phase  exists  in  each  reference 
control  volume  and  is  specified  solely  by  a  volume  fraction.  Such 
models  are  not  able  to  capture  the  limiting-current-density  phe¬ 
nomena. 

A  great  deal  of  the  recent  modelling  activity  has  been  concen¬ 
trated  on  multi-dimensional  numerical  simulation,  noteworthy 
examples  of  which  are  [2,4-7,9,18,19,27,29,30]  (by  no  means 
an  exhaustive  list).  The  multi-dimensionality  of  these  models 
significantly  complicates  their  solution,  forcing  approximations 
to  be  made  with  respect  to  water  and  the  CCL.  Natarajan  and 
Nguyen,  [5],  use  a  simplified  and  regularized  form  of  the  per¬ 
meability  in  their  isothermal  model  and  do  not  explicitly  model 
the  membrane  and  CCL.  He  et  ah,  [2],  adopt  a  similar  approach, 
simplifying  further  by  assuming  that  the  capillary  diffusion  co¬ 
efficient  is  constant.  Wang  et  al.,  [9],  describe  the  water  balance 
using  the  so-called  multi-phase  mixture  (M1 2 3 4)  model,  [31],  where 
thermodynamic  equilibrium  between  liquid  water  and  vapour  is 
assumed  to  prevail  throughout.  In  other  respects  their  model  is 
similar  to  those  in  [5,2],  An  isothermal  M2  model  is  also  used  in 
[30],  the  focus  of  which  is  flooding  in  the  GDL  and  CCL.  Siegel 
et  al.,  [6],  incorporate  the  CCL  and  membrane,  making  the  as¬ 
sumption  that  water  exists  only  as  vapour  or  as  a  species  dis¬ 
solved  in  the  electrolyte/membrane,  and  are  not  able  to  account 
for  liquid-water  effects  directly.  The  same  is  true  of  the  model 
in  [29],  which  treats  all  three  forms  of  water  as  a  single  phase. 
In  [7],  Siegel  et  al.  include  all  three  forms  but  introduce  several 
convoluted  assumptions  regarding  the  equilibrium  between  the 
three  phases,  as  well  as  a  highly  regularized  capillary  diffusion 
coefficient.  One-dimensional  models  that  rigorously  treat  water 
can  be  found  in  the  [3,32,33].  In  the  former,  Lin  et  al.  simplify 
the  capillary  diffusion  coefficient,  adopt  isothermal  conditions 
and  assume  a  constant  gas  pressure.  Nam  and  Kaviany,  [32], 
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concentrate  on  the  GDL,  treat  the  CCL  as  an  interface  and  only 
introduce  the  membrane  through  the  boundary  conditions.  They 
provide  a  thorough  and  rigorous  discussion  on  the  role  of  micro- 
structural  parameters  in  liquid-water  movement.  Pasaogullari  et 
al.,  [33],  employ  the  M2  model  in  a  detailed  study  of  water  ef¬ 
fects  in  the  membrane  and  GDL,  incorporating  both  anode  and 
cathode.  Their  model  is  isothermal  and  treats  the  catalyst  layer  as 
an  interface,  focussing  primarily  on  the  effect  of  a  micro-porous 
layer  attached  to  the  GDL. 

In  this  paper  we  develop  a  comprehensive  one-dimensional 
model,  explicitly  incorporating  a  GDL,  CCL  and  membrane. 
Most  of  the  features  separately  considered  in  detail  in  the  mod¬ 
els  described  above  are  included,  as  well  as  several  new  fea¬ 
tures,  such  as  a  separate  equation  for  dissolved  oxygen,  non¬ 
equilibrium  between  water  in  the  membrane  and  vapour,  and 
membrane  swelling.  Our  aim  in  this  paper  is  to  present  the  basic 
features  of  the  model,  and  to  demonstrate  its  success  in  captur¬ 
ing  the  convoluted  water-driven  phenomena.  Extension  of  the 
model  to  two  dimensions  and  to  one  of  the  entire  MEA  can  be 
achieved  in  a  straightforward  manner. 

The  outline  of  this  paper  is  as  follows.  In  the  next  section  we 
list  the  assumptions  adopted  in  deriving  the  model  and  discuss 
the  broader  issues  faced  in  MEA  modelling.  In  Section  3  we 
present  the  model,  with  complete  details  and  explanations.  In 
Section  4  we  describe  the  experimental  procedure.  The  design 
of  the  cell  is  such  that  variations  down  the  channel  are  minimal, 

1. e.,  pseudo  one-dimensional,  which  allows  comparison  between 
modelling  and  experimental  results  to  be  made.  The  numerical 
results  are  then  presented  in  Section  5,  in  which  we  first  demon¬ 
strate  agreement  with  the  experiments.  Finally,  in  Section  6  we 
provide  a  discussion  and  outline  future  work. 

2.  Model  description  and  assumptions 

Before  embarking  on  the  derivation  of  the  model,  it  is  in¬ 
structive  to  list  the  factors  that  are  commonly  believed  to  impact 
MEA  performance: 

•  the  diffusion  of  oxygen  through  the  pores  of  the  GDL  and 
CCL 

•  the  transport  of  oxygen  as  a  species  dissolved  in  liquid  water 
and  electrolyte 

•  water  production 

•  operating  temperature  and  temperature  variations 

•  the  micro-structure  of  the  GDL  and  CCL 

•  bulk  convective  motion  of  liquid  water  in  response  to  pressure 
gradients 

•  water  vapour  in  the  channels  and  in  the  pores  of  the  GDL  and 
CCL 

•  hydration  of  the  membrane 

•  electro-osmotic  drag  (migration) 

•  condensation/evaporation  and  absorption/desorption  (to  and 
from  the  electrolyte) 

•  the  transport  of  protons  and  electrons 

•  flooding  of  the  GDL  and  CCL  pores 

•  membrane  swelling 


•  non-uniform  proton  concentrations 

•  degradation  issues,  such  as  carbon  oxidation  and  hydrogen 
peroxide  production 

•  transient  effects. 

To  simplify  matters,  we  have  chosen  to  neglect  the  last  two 
items  in  this  work.  Their  incorporation  is  the  subject  of  current 
activity,  in  which  the  entire  MEA  is  explicitly  modelled.  We  say 
more  on  this  in  Section  6.  We  now  describe  the  model  and  list 
the  assumptions  adopted  in  its  derivation. 

1.  CCL  characterization.  As  in  the  class  of  models  discussed 
in  the  introduction,  we  assume  that  the  carbon  grains 
form  spherical  clusters  (agglomerates)  separated  by  primary 
pores.  Surrounding  the  agglomerates  is  an  electrolyte  layer, 
on  which  a  layer  of  liquid  water  can  also  exist.  We  assume 
that  the  contact  between  the  carbon  agglomerates  is  sufficient 
to  ensure  the  free  flow  of  protons  and  electrons,  and  that  the 
connecting  paths  posses  a  negligible  volume  fraction.  Within 
the  agglomerates  exists  a  fraction  of  the  electrolyte.  High- 
resolution  electron  microscopy  confirms  that  the  agglomer¬ 
ate  structures  are  coated  with  electrolyte  and  some  liquid 
water.  Middelman,  [34],  postulates  that  the  oxygen  does  not 
have  access  to  the  Pt  particles  located  within  the  agglom¬ 
erates,  particularly  those  close  to  their  centres.  For  reaction 
to  occur  within  the  agglomerates,  sufficient  electrolyte  must 
be  present,  the  quantity  of  which  depends  on  the  prepara¬ 
tion  process,  [10,12,35],  but  is  typically  low  (though  [  12,35] 
are  two  examples  of  efforts  directed  at  increased  electrolyte 
penetration  into  the  agglomerates).  We  therefore  assume  that 
the  Pt  within  the  agglomerates  is  not  active,  i.e.,  that  it  does 
not  contribute  to  the  bulk  reaction.  It  is  possible  to  derive 
an  effectiveness  factor  to  remove  this  assumption.  However, 
since  it  is  unclear  that  significant  reaction  occurs  within  the 
agglomerates,  and  for  the  sake  of  simplicity,  we  continue 
with  our  assumption.  If  needed  or  justified,  the  inclusion  of 
an  effectiveness  factor  is  straightforward. 

2.  Oxygen  transport.  Diffusion  in  large  pores  is  assumed  to 
be  predominantly  continuum,  given  their  characteristic  size 
and  the  typical  operating  pressures.  Assuming  that  the  oxy- 
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gen  and  water  vapour  are  small  concentrations  in  a  nitrogen 
gas,  we  employ  Fick’s  law.  The  oxygen  in  the  primary  pores 
reaches  the  Platinum  particles  on  the  agglomerate  surface  by 
dissolving  in  the  water/electrolyte  and  diffusing  across  the 
layer(s).  Henry’s  law  is  used  to  describe  the  equilibrium  in 
oxygen  concentration  at  the  interface  between  the  gas  and 
liquid/solid  phases  (see  Fig.  2). 

3.  Non-uniformity  of  proton  concentration.  The  protons  are 
transported  in  the  form  of  hydronium  ions,  mO-1- .  We  assume 
electro-neutrality  (the  number  of  negative  charges  equals  the 
number  of  positive  charges)  but  we  do  not  make  the  typical 
assumption  of  a  homogeneous  distribution  of  charged  sites 
in  the  electrolyte,  which  would  allow  us  to  take  the  concen¬ 
tration  of  H30+  as  constant. 

4.  Water.  The  model  accounts  for  water  in  all  three  forms;  as 
a  dissolved  species,  as  vapour  and  as  liquid.  Liquid  water 
resides  in  the  primary  pores.  The  water  dissolved  in  the  elec¬ 
trolyte  is  transported  by  diffusion  (concentration  gradients) 
and  electro-osmotic  drag  (migration).  As  do  Seigel  et  al., 
[6],  we  assume  that  the  net  water  produced  is  either  in  liq¬ 
uid  or  vapour  form.  Condensation  and  evaporation  are  mod¬ 
elled  using  the  approach  in  [2,4,5],  and  references  therein. 
In  this  representation  the  transfer  occurs  at  a  finite  rate,  dic¬ 
tated  by  the  deviation  of  the  local  thermodynamic  state  from 
equilibrium.  In  a  similar  fashion,  we  introduce  phase  change 
between  vapour  and  dissolved  water  by  considering  the  devi¬ 
ation  from  equilibrium  between  dissolved  water  and  vapour. 
Details  for  each  are  presented  later. 

5.  Electrolyte  volume.  We  include  the  effects  of  swelling  on  the 
volume  of  electrolyte  in  the  CCL. 

6.  Temperature.  The  model  is  non-isothermal,  which  is  neces¬ 
sary  to  model  phase  change  accurately  and  to  study  the  effects 
of  water  activity.  We  do  however  assume  a  single  temperature 
for  all  phases. 

Other  assumptions  are  (a)  negligibly  small  viscous  flows, 
[36],  and  (b)  a  negligible  ohmic  potential  drop  in  the  electron¬ 
ically  conducting  carbon  phase,  justified  by  the  relatively  high 
carbon  conductivity. 


nine,  T,  CD  +hpe(HCp  -  ce). 


reactant  consumption 


(2) 


d 

dx 


((1  -  s)Ds 


=  0, 


(3) 


where  cc  is  the  volume  fraction  of  electrolyte  and  Dp,  Z>n  and  Dt 
are  the  molecular  diffusion  coefficients  of  oxygen  and  nitrogen 
through  the  primary  pores  and  of  oxygen  through  the  electrolyte. 
The  first  are  subject  to  a  Bruggeman  correction  (in  the  absence  of 
a  value  for  the  tortuosity  of  the  flow  path)  and  their  dependence 
on  temperature  and  pressure  follows  the  Chapman-Enskog  for¬ 
mula  (in  m2  s_1),  [37], 


7^3/2  7^3/2 

Dp  =  3.8  x  1(T9— e]!2,  DN  =  3.6  x  1(T9—  e3/2,  (4) 

where  ep  is  the  volume  fraction  of  primary  pores,  T  is  temper¬ 
ature  and  P  is  the  gas  pressure.  The  variations  of  Dp  and  £>n 
with  temperature  and  pressure  are  neglected,  that  is,  we  use  a 
constant  value  based  on  the  temperature  at  the  cathode  channel. 
Since  temperature  and  pressure  variations  are  weak,  the  corre¬ 
spondingly  small  changes  in  Dp  and  have  a  negligible  effect. 
The  coefficient  De  has  the  following  form  (taken  form  [38]): 

De  =  3.1  x  10“7e“2768/:r,  (5) 


where  once  again  we  approximate  temperature  in  this  expression 
by  its  value  at  the  cathode  channel. 

7l(r]c,  T,  C*)  is  the  reaction  rate,  based  on  the  Butler- Volmer 
law,  /ipe  is  a  volumetric  mass  transfer  coefficient  from  gas  to 
electrolyte  (on  the  air  side),  and  H  is  a  dimensionless  Henry 
constant.  The  diffusive  flux  of  oxygen  from  the  pore  space  to 
the  interface  between  the  gas  and  electrolyte/water  is  balanced 
by  the  amount  adsorbed  into  the  electrolyte/water. 

To  prescribe  the  mass-transfer  coefficient,  h  pe ,  we  use  the 
formula  for  the  Sherwood  number  given  in  [37]  for  flow  past  a 
spherical  particle. 

Sh  =  2  +  0.6  Rel/3Sc2/3, 


3.  Model  derivation 

3.1.  Catalyst  layer 

The  geometry  is  depicted  in  Fig.  1,  showing  x  as  the 
space  variable,  with  the  cathode  channel  placed  at  x  =  0,  the 
GDL/CCL  interface  at  x  =  L  i ,  the  membrane/CCL  interface  at 
x  —  L2  and  the  membrane/ ACL  interface  at  x  =  L3. 

Let  5  be  the  saturation  (fraction  of  primary  pore  volume  occu¬ 
pied  by  liquid  water),  the  concentrations  of  oxygen  and  nitrogen 
in  the  primary  pores  be  Cp  and  Cn,  respectively,  and  the  concen¬ 
tration  of  oxygen  in  the  electrolyte  be  Ce.  Taking  into  account 
reaction,  diffusion  and  absorption,  mass  balances  yield: 

d  (  dCp\ 

( (1  -  -vjDp  J)  -  -  hpfllCp  -  Ce)  .  (1) 


where  Re  is  the  Reynolds  number  and  Sc  is  the  Schmidt  num¬ 
ber.  For  a  predominantly  diffusive  process,  we  can  approximate 
Sh  =  2.  From  the  formula  hpesa  =  ShDp/Dag,  where  Dd„  is  the 
agglomerate  diameter  and  sa  is  the  specific  surface  area  of  the 
agglomerates,  we  obtain  hpe  =  (?(105).  It  should  be  noted  that 
estimating  mass  transfer  coefficients  is  itself  a  great  challenge,  a 
detailed  discussion  of  which  is  beyond  the  scope  of  the  present 
paper. 

The  electrolyte  volume  fraction,  ee,  will  increase  as  it  swells 
with  water.  We  separate  the  electrolyte  film  that  coats  the  ag¬ 
glomerates  (volume  fraction  ejj  from  that  contained  in  the  inte¬ 
riors  (volume  fraction  c'c)  and  write  e9  =  The  Pt  inside 

the  agglomerates  is  assumed  to  be  inactive,  that  is,  there  is  a 
negligible  surface  area  of  Pt  in  contact  with  the  fraction  c'c  of 
electrolyte.  The  combined  volume  fraction  of  the  carbon,  Pt  and 
small  pores,  ea,  is  assumed  constant.  The  volume  fraction  of  pri¬ 
mary  pores  is  then  given  by  ep  =  1  —  ea  —  ee.  The  swelling  of 
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the  electrolyte  is  assumed  to  impact  on  the  thickness  of  the  elec¬ 
trolyte  film.  To  incorporate  it  we  use  the  following  approximate 
relationship: 

ee  =  el  +  ksX.,  (6) 


where  ks  =  0.0126,  and  X  is  the  membrane  water  content  (mol 
H20/mol  SO^~),  given  by 


1  —  ksC(\ 


(7) 


where  C d  is  the  dissolved  water  concentration  in  the  electrolyte, 
normalized  with  respect  to  v  —  1800  mol  m-3,  the  fixed-charge 
site  concentration  of  the  membrane,  [39]. 

To  derive  an  equation  for  the  electrolyte  potential  in  the  CCL, 
(j>,  we  consider  a  mass  balance  for  H^O3-.  With  the  assumption 
of  a  negligible  convective  flux,  the  contribution  to  the  movement 
comes  from  migration  and  concentration  gradients: 


d 

dx 


eecre  d(f>  dCH\ 

F  dx  fe  "  d  dx  ) 


diffusive  flux  of 


K(r,c,  7]  Cl)  , 

v. _ v _ ✓ 

H3  0+ consumption 


(8) 


where  C h  is  the  FFO-1-  concentration,  cre  is  the  protonic  conduc¬ 
tivity  of  the  electrolyte,  D\  \  is  the  H30+  diffusion  coefficient  and 
F  is  Faraday’s  constant.  Ch  and  Cd  are  related  as  follows  [40] 


1 

Ch  =  —~k&C, d  + 


+  keCd, 


where  the  temperature  dependent  ke  is  given  by 
k'  =  t  °exp(? 


(9) 

(10) 


As  an  approximation  we  assume  that  kt  is  constant.  Note  that 
Eq.  (9)  is  derived  with  the  assumption  of  electro-neutrality.  We 
refer  to  [41]  for  details.  The  diffusion  coefficient  D\\  is  taken 
from  [40]: 

Dh  =  1.6  x  10“8e“1683/:rCd,  (11) 


where  once  again  we  neglect  the  dependence  on  temperature 
variations.  For  ere  we  use  the  law  derived  in  [39], 

ae  =  exp  ^1268  ^  (5.14A.  -  0.326).  (12) 

The  algebraic  term  on  the  right-hand  side  of  (12)  accounts 
for  the  effect  of  hydration  on  the  conductivity.  It  is  based  on 
empirical  data  collected  by  Springer  et  al.  [39]. 

An  equation  for  the  temperature,  T,  is  derived  from  an  energy 
balance  of  conduction,  convective  heat  flux  from  liquid  water 
movement  and  heat  sources 

d  /  dT  \  /  (—8s)T  \ 

-  —  Ik—  -  sevp\C\V\T\  =  ( — - - h  Fr/J  T,  Ce) 

v - v - ' 

reaction  and  activation  loss 

+  hg\hvc{RTCy  -  P Sat) .  (13) 

" - v - " 

phase  change 

ohmic  heating 


'  d<p 

+  w  1  s 


evaporation  condensation 

RH<1  /  \  RH>1 


Fig.  3.  A  schematic  representation  of  the  interactions  between  water  in  the  three 
phases. 

In  this  equation,  —  (8s)  is  the  entropy  associated  with  ORR, 
C j  is  the  water  vapour  concentration,  ne  —  4  is  the  number  of 
electrons  transferred  in  one  reaction,  /jpc  is  a  mass  transfer  co¬ 
efficient,  defined  below,  h„}  is  the  liquid-gas  enthalpy  change,  p\ 
is  the  density  of  liquid  water,  C\  is  the  specific  heat  capacity  of 
liquid  water,  and  k  is  the  effective  thermal  conductivity,  found 
from  a  volume  average  of  the  individual  conductivities.  For  sim¬ 
plicity,  we  use  constant  values  for  k.  Finally,  the  quantity  i>i  is 
the  liquid-water  interstitial  velocity,  which  we  define  later. 


3.1.1.  Water  in  the  catalyst  layer 

The  schematic  in  Fig.  3  depicts  the  features  of  the  liquid- 
water  model  and  the  assumptions  we  make  in  deriving  it.  For  a 
detailed  description  of  liquid-water  transport  in  MEA  we  refer 
to  Nam  and  Kaviany  [32]. 

The  mass  balance  for  water  dissolved  in  the  electrolyte  is 
written  as  follows: 


d 

dx 


/ 


dQ 

dx 


V 


5  d0 

44  Fv  dx 


drag 


-  /!dv(Cd  -  C*), 

■s - v _ ✓ 

vapour-o-dissolved 


(14) 


where  Dd  is  the  diffusion  coefficient  for  dissolved  water.  The 
water  vapour  mass  balance  is  expressed  as  follows: 

d  (  dcr\ 

-  —  (^(1  -  s)Dv  J  J  =  -  hpc(RTCy  -  Rsat) 

liquids  gas 

+  hdvv(Cd  -  c*d)  +  y ftfac,  7]  Cl),  (15) 

^  -  s 

H2O  production 


where  Dv  is  the  vapour  diffusion  coefficient.  The  switch  function 
r  satisfies: 


r  = 


0  a^j  1 
1  v  tX  1 


(16) 


The  bulk  diffusivity  Dv  is  obtained  from  Chapman-Enskog  the¬ 
ory  and  a  Bruggeman  correction,  and  Dd  from  the  relationship 
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provided  in  [42]  (with  temperature  and  pressure  approximated 
by  their  values  at  the  cathode  channel): 


Dd  =  3.5  x  l(T5ke“2436/:r, 


Dv=  4.1  x  10“9 


(17) 


For  liquid  water,  the  mass  balance  is 

-  ^  (^m)  =-V(flrcv-psat)+  t,  d), 

"capillary  action"  H2°  production 

(18) 


where  i>i  and  W\  are  respectively  the  interstitial  velocity  and 
molar  mass  of  liquid  water.  v\  is  determined  by  the  Darcy-law 
approximation  to  momentum  conservation: 


ki  dp 


v\  =  —  —  , 
pi  1  d.r 


(19) 


where  k\,  /i\  and  p  are  respectively  the  relative  permeability, 
viscosity  and  pressure  of  the  liquid.  The  latter  is  given  by  p  = 
P  —  pc,  so  that  differentiating  Eq.  (18)  and  using  Eq.  (19)  yields: 


Pi  d  f  I"  dpc  ds 

d.v  V1’'0  5  [  ds  dx 


d  P 

d.v 


— —hpC(RTCv  —  /"’sat) 


3(r  —  i) 

+  2  JTZ(Vc,T,Cl).  (20) 

There  is  a  great  deal  of  debate  surrounding  both  the  form  and 
magnitude  of  the  permeability  and  capillary  pressure.  Mazumder 
and  Cole  [4],  compare  the  Leverett  model  employed  in  [5]  with 
the  empirical  relationship  in  [9] ,  in  which  both  of  these  quantities 
assume  entirely  different  forms,  yielding  markedly  different  pre¬ 
dicted  levels  of  saturation.  There  are  also  other  models  of  liquid 
transport  in  porous  media,  of  varying  degrees  of  sophistication. 
For  example,  Nam  and  Kaviany,  [32],  scale  saturation  against  an 
immobile  value  (below  which  the  liquid  water  is  discontinuous 
and  therefore  does  not  move),  and  place  the  scaled  saturation 
inside  the  Leverett  function.  All  of  the  forms  mentioned  can  be 
easily  implemented  in  our  model  but  we  choose  to  follow  that 
in  [4,30,33]  (and  references  therein). 

Pc  =  crcJ(s),  k\(s)  =  kcs3,  (21) 


where  kc  is  the  absolute  permeability  of  the  CCL  and  J(s)  is 
the  capillary-pressure  function.  The  quantity  crc  is  defined  as 
follows: 


where  0LC  is  the  contact  angle  for  the  CCL  (0  <  0LC  <  n/2)  and 
a'  is  the  surface  tension.  The  grouping  2y/ep//rc  is  then  equal 
to  the  characteristic  capillary  radius. 

The  Leverett  function  for  a  hydrophilic  medium,  for  which 
the  wetting  phase  is  the  liquid,  is  given  by 

J(s)  =  1.417(1  -  s)  -  2.12(1  -  sf  +  1.262(1  -  sf.  (23) 


The  liquid-vapour  phase-change  term,  derived  from  that 
given  in  [4]  and  references  therein,  is  driven  by  the  deviation 


from  equilibrium,  RTCV  —  Ps at,  where  PsM  is  the  saturation  pres¬ 
sure  of  water  and  the  first  term  is  equivalent  to  the  partial  pressure 
of  the  vapour.  The  latter  is  defined  by  XVP ,  where  P  is  the  pres¬ 
sure  of  the  gas  phase  and  Xv  is  the  molar  fraction  of  vapour. 
Using  the  ideal  gas  law,  P  =  CRT,  XY  is  given  by 


Cv  CVRT 
~~C  ~  P 


(24) 


where  C  —  Cv  +  Cn  +  Cp  is  the  molar  density  of  the  gas.  A 
relationship  for  the  saturation  pressure  is  provided  in  the  discus¬ 
sion  of  the  boundary  conditions.  In  order  to  distinguish  between 
evaporation  and  condensation,  the  mass-transfer  coefficient  hpc 
depends  on  the  sign  of  the  driving  force  R  TCV  —  PsM  and  its 
precise  form  is 


/zpc  — 


Kc€p(lRTXv  RTCV  >  Psat 
K/ff  RTCV  <  Psat  ’ 


(25) 


where  atc  and  zcc  are  the  condensation  and  evaporation  rate  con¬ 
stants.  The  values  of  ke  and  kc  are  taken  from  [2],  We  use  a  con¬ 
tinuously  differentiable  implementation  of  the  following  form: 


^  Kc£p(l  —  ■'>)XV  /  I R TC\  —  .Psatl  \ 

'pc  “  2RT  V  +  RTC,  -  Lsat  ) 

i  zceep.tpi  /  I RTCW  —  Fsat|  \ 

2 Wi  V  fiTCv  -  Psat  )  ‘ 

In  a  similar  fashion,  the  vapour-dissolved  phase-change  term 
inEqs.  (14)  and  (15)  is  driven  by  the  deviation  from  equilibrium, 
Cd  —  C^,  where  C(j  is  the  concentration  of  dissolved  water  in 
equilibrium  with  vapour.  The  latter  is  derived  from  the  following 
law  for  the  equilibrium  water  content  measured  at  80  °C  (found 
in[43]): 


k*  =  0.3  +  10.8tzw  -  16 al  +  U.lal, 

X* 

c;  =  - ,  (27) 

d  1+  0.0126k* 

where  czw  is  the  water  vapour  activity.  To  distinguish  between 
desorption  and  water-uptake  of  the  electrolyte  (absorption),  the 
mass-transfer  coefficient  hdv  depends  on  the  sign  of  the  driving 
force  Cd  —  Cj .  From  the  experimental  and  numerical  results 
in  [44],  we  can  approximate  the  coefficients  of  absorption  and 
desorption  as  follows: 


/*dv  — 


Ka(l  -  S)X  Cd  -  C*  <  0 
*d(l  -  S)X  Cd  -  C*  >  0  ’ 


(28) 


where  zq  and  zce  are  given  in  Table  1 .  Their  weak  dependence 
on  temperature  variations  is  neglected.  The  factor  1  —  s  ac¬ 
counts  for  the  blockage  of  the  pores  by  the  liquid  water.  To 
implement  this  approach  we  use  a  continuously  differentiable 
form  of 

=  *„<1  -  S)X  (l  +  ]^^1) 

<29) 


Table  1 
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Symbol 

Quantity 

Size 

References 

L 

Catalyst  layer  thickness 

22.5  |xm 

L  m 

Membrane  thickness 

62.5  |xm 

Lg 

Cathode  GDL  thickness 

200  p,m 

L  (Tc) 

Anode  (cathode)  channel  temperature 

70  °C 

av/,c  (flw,a) 

Cathode  (anode)  channel  water  activity 

1(1) 

c 

Total  gas  concentration  in  cathode  channel 

Pc/(RTc)mo\  m-3 

Cp 

Oxygen  concentration  in  cathode  channel 

0.05C  mol  m-3 

Pc 

Gas  pressure  in  the  cathode  channel 

30  psig  (206.842  kPa) 

6P 

Volume  fraction  of  primary  pores 

0.6 

[10-12] 

4 

Electrolyte  volume  fraction  in  agglomerates 

0.2 

[10-12] 

Volume  fraction  of  C,  Pt  and  small  pores 

0.3 

[10-12] 

6g 

Porosity  of  the  GDL 

0.75 

[49] 

CO2,  ref 

Reference  oxygen  concentration 

4.55  mol  m~3 

[6] 

A 

Oxygen  diffusivity  in  liquid  water  (60  °C) 

4.82  x  10-9  m2  s-1 

[37] 

<*C 

Cathodic  transfer  coefficient 

0.55 

Anodic  transfer  coefficient 

0.45 

E0 

Open  circuit  potential 

0.95  V 

A g 

Agglomerate  diameter=  2/?ag 

0.4  p,m 

[10-12] 

<5o 

Electrolyte  film  thickness  without  swelling 

15  nm 

[17] 

N 

Number  of  agglomerates  per  unit  volume 

1.257  x  1019  m~3 

Est. 

apt 

Specific  surface  area  of  Pt 

1000  cm2  (mg  Pt)— 1 

[50] 

mp  t 

Pt  loading 

0.4  (mg  Pt)  cm~2 

[3] 

hpe 

Oxygen  mass  transfer  rate 

10s s"1 

Est. 

H 

O2  Henry’s  law  constant  (dimensionless) 

0.3 

[51] 

Ka  (<Cd) 

Absorption  (desorption)  constant 

10_5/9  (10_5/3)  ms-1 

[2] 

4C  (4S) 

CCL  (GDL)  contact  angle 

90°  (120°) 

([33]) 

Kc  (Kg) 

Absolute  permeability  of  CCL  (GDL) 

10~13  (8.7  x  10-12)  m2 

[52]  ([49]) 

Ml 

Liquid- water  viscosity 

10~3  KgnG1  s_I 

[4] 

o' 

Surface  tension 

0.07  Nm”1 

[4] 

k 

Thermal  conductivity  of  CCL 

2Wm_1  Kr1 

[29] 

km 

Thermal  conductivity  of  membrane 

lWm-'r1 

[29] 

ka 

Thermal  conductivity  of  GDL 

4Wm-‘r1 

[29] 

Liquid- water  removal  constant 

0.075  itG1 

Fitted 

102,c 

Reference  cathode  exchange  current  density 

40  A  m“2 

Fitted 

PiC 1 

Heat  capacitance  of  water 

4.184  x  106  J  m-3  K-1 

Est. 

-8s 

Entropy  associated  with  ORR 

162.2  Jmor1  K_1 

[53] 

3.1.2.  Schroeder’s  paradox 

It  is  well  known  that  the  water  content  X  depends  on  the  water 
activity,  aw  =  Xv  P/  PsM.  The  precise  relationship  was  correlated 
by  Hinatsu  et  al.  in  [43], 

kv  =  ka  =  0.3  +  I0.8flw  -16 al+  14.  la3,  (at  353  K),  (30) 

which  has  a  maximum  of  Xv  =  9.1  at  equilibrium  with  vapour. 
However,  when  the  electrolyte  is  submerged  in  liquid  water  its 
equilibrium  water  content  appears  to  jump  discontinuously  to  a 
higher  value,  which,  following  [44],  we  take  to  be  X  —  16.8. 

In  an  attempt  to  capture  this  anomaly,  known  as  Schroeder’s 
paradox,  we  define  a  critical  saturation,  s*,  at  which  the  agglom¬ 
erates  are  entirely  coated  with  liquid  water  and  equilibrium  with 
liquid  water  is  achieved.  This  critical  value  is  approximated  as 
the  immobile  saturation,  below  which  the  water  mass  is  not  con¬ 
tinuous  (and  therefore  does  not  flow).  We  point  out  that  this  is 
an  approximation,  based  on  the  assumption  that  the  electrolyte 
is  continuously  covered  by  water  for  s  >  s*.  We  use  the  value 
of  s*  —  0.1  given  in  [32],  and  refer  to  that  paper  for  references 
related  to  its  measurement.  The  water  content  at  the  boundary 
between  the  CCL  and  membrane  is  kept  continuous  by  enforc¬ 


ing  equilibrium  between  the  membrane  and  liquid  water  when 
s(L,2)  =  s*.  Alternatively,  we  may  assume  that  the  water  content 
in  the  membrane  drops  continuously  to  a  vapour-equilibrated 
value  at  x  —  L-\,  but  this  does  not  qualitatively  alter  the  results. 
We  therefore  make  the  simpler  choice  and  avoid  further  assump¬ 
tions  in  specifying  the  rate  of  decrease. 

A  list  of  parameters  and  their  values  can  be  found  in  Table  1. 


3.2.  Reaction  rate  and  limiting  current  density 


The  ORR  rate  (in  mol  m  3  s  1 )  is  given  by  the  Butler- Volmer 
law  and  first-order  kinetics  in  oxygen  concentration: 


K(rlc,T,Cl)  = 


al02.c  r 
tp  ee°e 

f'c02.ref 


oiiFric/RT 


e-acFVc/RT  j  ;  (31) 


where  i q2  c  is  the  cathode  exchange  current  density,  aa  andac  are 
the  anodic  and  cathodic  transfer  coefficients,  co,  ,ref  is  a  reference 
oxygen  molar  concentration,  a  is  the  surface  area  of  catalyst  per 
unit  volume  of  catalyst  layer,  r/c  is  the  overpotential  and  R  is 
the  universal  gas  constant  (8.314  Jmol-1  K-1)).  The  quantity  a 
is  a  function  of  the  specific  Pt  surface  area  (Pt  surface  area  per 
unit  mass  of  Platinum),  apt,  the  Pt  loading,  mpt,  and  the  CCL 
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thickness,  L: 
a  =  cipttnpt/L. 

However,  the  active  platinum  surface  area  is  much  more  diffi¬ 
cult  to  characterize.  The  overpotential,  r]c,  is  defined  precisely 
through  the  relationship: 


r]c  =  Uc  -  E0  -  </>, 


(32) 


where  Uc  is  the  potential  of  the  carbon  matrix  (constant  by  as¬ 
sumption)  and  Eq  is  the  so-called  open  circuit  potential,  as¬ 
sumed  constant.  The  value  of  oxygen  concentration  in  the  Butler- 
Volmer  expression,  (31),  ought  to  reflect  accurately  that  reaction 
takes  place  at  the  surfaces  of  the  carbon  agglomerates.  The  con¬ 
centration  at  these  surfaces,  C*  ,  is  generally  different  from  the 
bulk  value  in  the  electrolyte,  particularly  given  the  restricted 
diffusion  of  oxygen.  C*  can  be  related  to  the  bulk  value,  Ce, 
by  balancing  the  rate  of  reaction  with  the  rate  of  diffusion  of 
oxygen  to  the  surface  of  the  agglomerates  (at  steady  state).  This 
mass  balance  can  be  approximated  as  follows: 

y'(Ce  -  Cl)  =  l-lZ{r)c,  T,  C“),  (33) 

where  y'  is  the  rate  of  oxygen  diffusion  through  the  elec¬ 
trolyte/water  film  to  the  surface  of  the  agglomerates  (in  s-1). 
Using  the  definition  of  T,  Csc)  and  solving  the  resulting 

equation  for  C*  then  yields: 


Cs  =  _ /Ce _ 

e  yr  ~ \~  €&/  (oPaFtlc/ RT  —  Q—McFric/ RT^J 

^02,  C 


where 


4Tco2,ref 


(34) 


so  that  the  final  form  of  the  reaction  rate  is 


nine,  t,  cl)  =  nine,  t,  ce) 

y'eerCe  (e?*F,ic/RT  -  q-^f^/rt^ 

=  /  +  eer  (e“a^c/«7'  _  e-«T'ic /*r) '  (35) 

In  order  to  relate  the  diffusion  rate  y'  to  the  microscopic 
properties  of  the  CCL,  we  define  an  electrolyte  film  thickness, 
<5e,  an  agglomerate  radius,  Rd g,  and  the  number  of  agglomerates 
per  unit  volume,  N: 

.  ,  De  ,  ADe 

y  (Ce  —  CD  &  A  -f(Ce-CD  so  that  Z=-r1,  (36) 

v - v _ ✓ 

molar  flux 


The  electrolyte  swelling  results  in  a  volume  change  equal  to  (per 
agglomerate)  (?[  —  efy/N.  The  electrolyte  film  thickness  is  then 
given  by 

&e  =  (^ag  +  <^)3  +  3(4JrV°))  -  RV  (38> 

When  s  >  s*,  as  previously  defined,  the  liquid  water  present 
in  the  catalyst  layer  will  provide  additional  resistance  to  the 
oxygen,  with  a  different  diffusion  coefficient  D\.  Given  the  hy¬ 
drophilic  nature  of  the  electrolyte,  we  assume  that  any  liquid 
water  coats  the  entire  surface  of  the  agglomerates.  Based  on 
these  assumptions,  the  thickness  of  the  water  layer  is  given  by 

S1=  ^ag  +  4)3  +  -^j  -i?ag,  (39) 

with  a  total  film  thickness 


S  =  8e  +  5i 


(40) 


The  quantity  y'  then  has  to  be  modified;  it  is  decomposed  into 
a  term  arising  from  diffusion  through  the  water  layer,  y[,  and  a 
term  of  the  type  described  above,  y'e.  To  approximate  the  concen¬ 
tration  of  oxygen  at  the  agglomerate  surfaces  two  balances  need 
to  be  performed,  one  for  diffusion  through  the  water,  to  relate 
the  bulk  concentration  to  concentration  at  the  water/electrolyte 
interface,  and  one  between  diffusion  through  the  electrolyte  and 
reaction,  to  relate  the  latter  concentration  to  C*.  The  reaction 
rate  then  takes  the  form  (35)  with 


II 

+  g" 

CD  ■>. 

,  ,  A’Di 

where  =  - , 

<5i 

AD„ 

) 

<D 

co 

II 

A'  =  4it{Rag  +  Se)2N. 

The  Henry  constant  H  would  likewise  require  modification.  It 
is  approximately  0.3  for  water  under  typical  operating  condi¬ 
tions,  and  was  assigned  the  value  0.15  for  Nafion  in  [27].  For 
simplicity,  we  assume  the  water-based  value  throughout  and  ap¬ 
proximate  it  as  a  constant. 

Note  that  the  quantity  IV  is  a  function  of  both  the  radius  of 
the  agglomerates  and  their  distribution,  for  example,  the  more 
densely  packed  the  agglomerates  the  larger  N. 

The  current  flow  is  the  reverse  of  the  H30+  flow,  so  that  the 
current  density,  I'  Am-2,  is 


y'eeCe  (ea'Fr>c/RT  -  e-a^F^RT) 
- - - - —  dx 

y'  -)-  ^QaaFr]c/RT  —  Q—acFr]c/ RT'j 


(42) 


where  A  is  the  specific  surface  area  of  agglomerates  (we  assume 
that  all  of  the  surface  area  is  covered,  whereupon  A  =  4i rR^N). 
The  volume  of  electrolyte  attached  to  the  surface  of  each  agglom¬ 
erate  is  el/N,  which,  from  our  assumptions,  covers  the  entire 
surface  of  the  agglomerate.  The  thickness  (without  swelling) 
and  6q  are  therefore  related  as  follows: 


S 


e  _ 

o  — 


4: rN  J 


t/3 


R 


ag- 


(37) 


3.3.  Membrane 


We  now  describe  the  equations  employed  in  the  membrane, 
L2  <  x  <  L3,  within  which  we  assume  that  oxygen  and  vapour 
do  not  exist,  i.e.,  the  membrane  is  impenetrable  to  both.  The 
membrane  potential  is  modelled  in  a  manner  similar  to  the  elec¬ 
trolyte  potential  in  the  CCL,  Eq.  (8), 


d 

dx 


/  cre  d 0 
\  F  dx 


+  -DftCd 


=  0, 


(43) 
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where  ere  and  £>h  were  previously  defined.  The  movement  of 
water  in  the  membrane  is  either  as  a  dissolved  species  or  pos¬ 
sibly  by  convection  as  a  liquid,  the  latter  depending  on  the 
state  of  the  membrane.  It  has  been  postulated  by  Weber  and 
Newman,  [45],  that  the  convective  movement  increases  with  in¬ 
creased  water  content  because  the  membrane  walls  are  inflated 
by  the  action  of  the  water.  However,  as  Mazumder  demon¬ 
strates  in  [46],  the  role  of  convection  is  typically  outweighed 
by  both  diffusion  and  electro-osmotic  drag,  even  when  a  pres¬ 
sure  increase  of  2  atm  from  anode  to  cathode  is  maintained.  We 
therefore  neglect  it,  particularly  since  we  will  not  consider  such 
substantial  pressure  drops.  The  equation  for  dissolved  water  is 
then 


5  d  0\ 

- Acre  — 

44  Fv  dx  ) 


=  0. 


(44) 


Finally,  the  equation  for  liquid  water  is 


Pi  d 
fi\W\  dx 


eg/ciCs) 


d  pc  ds 
ds  dx 


d  P 

dx 


— —  hpC(RTCv— Tsat), 

(49) 


with 


pc  =  ogJ{s),  K\(s)  =  Kgs3,  erg  =  a'  cos0f,  — , 

Kq 


(50) 


where  Kg  is  the  absolute  permeability  of  the  GDL  and  0§  is 
the  contact  angle  (jr/2  <  0§  <  it).  The  Leverett  function  for  a 
hydrophobic  medium,  for  which  the  wetting  phase  is  the  gas 
phase,  is  given  by 


J(s)  =  1.4175  -  2.12s2  +  1.262s3. 


(51) 


The  equation  for  temperature  now  involves  only  a  conduction  3.5.  Boundary  conditions 
term  and  ohmic  losses 


d  f  d T\ 
dx  \  m  dx  / 


(45) 


where  km  is  the  effective  thermal  conductivity  of  the  membrane. 
At  the  interface  between  the  membrane  and  CCL  a  further  con¬ 
dition  is  enforced  to  ensure  the  continuity  of  the  water  content, 
X.  This  condition  was  discussed  above. 


3.4.  Gas  diffusion  layer 


In  the  GDL,  0  <  x  <  L  1 ,  we  specify  balances  for  the  masses 
of  oxygen,  nitrogen,  vapour  and  liquid  water,  together  with  an 
energy  balance.  Equations  for  the  first  two  are 


d 

dx 


d-i) 


d 

dx 


Z?g(l-s)^ 

dx 


=  0, 


(46) 


where  Dp  and  i/  are  the  effective  molecular  diffusion  coeffi¬ 
cients  for  oxygen  and  nitrogen  in  the  GDL,  and  eg  is  the  GDL 
porosity. 

The  equation  for  vapour  concentration  is 
d  (  rdCT\ 

f  (1  -  s)D^j  =  -hvc(RTCY  -  P sat), 


For  the  discussion  of  the  boundary  conditions  we  recall  Fig.  1 , 
which  depicts  the  geometry.  At  the  interface  between  the  mem¬ 
brane  and  CCL,  x  =  Li,  the  fluxes  of  oxygen  (in  the  pores), 
water  vapour,  nitrogen  and  liquid  water  are  taken  to  be  zero; 
that  is,  these  species  are  assumed  not  to  penetrate  the  mem¬ 
brane.  It  is  possible  that  small  amounts  of  the  dissolved  oxygen 
reach  the  anode  and  small  amounts  of  hydrogen  reach  the  cath¬ 
ode.  These  are  primarily  degradation  issues,  which  are  dealt 
with  in  a  forthcoming  paper.  Since  protons  have  no  effective 
transport  mechanism  in  the  GDL,  the  flux  of  protons  at  the  in¬ 
terface  between  the  GDL  and  CCL,  x  =  L 1,  is  zero.  Similarly, 
the  dissolved-water  and  oxygen  fluxes  must  be  zero  at  x  =  L 1 
in  order  to  ensure  their  continuity: 


dCe 

dx 


_  ^  dCd  ,  5  '  d (p_d(p 

d  dx  44  Fv  e  dx  dx 
dCp  dCe  dCv  ds 

dx  dx  dx  dx 


=  0,  x—L  1 
x  =  L2. 


(52) 


At  the  interface  between  the  cathode  channel  and  the  GDL 
interface,  x  =  0,  the  oxygen  and  vapour  concentrations  are  pre¬ 
scribed.  At  the  membrane/ ACL  interface,  x  =  L  3,  the  membrane 
potential  is  approximated  as  zero 


Cp(0)=Cp,  Cv(0)  =  Cv,  <KLi)  =  0.  (53) 


(47) 

where  l)\f  is  the  effective  molecular  diffusion  coefficient  for 
vapour  in  the  GDL.  Df  l/  and  D*)  are  approximated  as  in 
(4)  and  (17)  (with  ep  replaced  with  eg),  evaluated  at  the  cathode 
channel  conditions.  The  phase-change  term  on  the  right-hand 
side  of  (47)  is  defined  as  in  (25)  and  (26),  with  ep  replaced  with 

eg- 

The  energy  equation  balances  conduction,  convective  heat 
flow  arising  from  the  liquid  water  motion  and  heating  from  phase 
change. 


The  value  of  Cv  is  calculated  from  the  water  activity  of  the 
cathode  channel,  awc,  which  is  defined  as  XYPC/Psa t,  where 
Zv  is  the  molar  fraction  of  vapour  in  the  cathode  channel,  Pc 
is  the  cathode-channel  gas  pressure  and  Psm,c  is  the  cathode- 
channel  saturation  pressure.  The  saturation  pressure,  a  function 
of  temperature,  is  given  by  the  formula  in  [39]: 

log10  Psat  =  -2.1794  +  0.02953(7’  -  273)  -  9.1837 

x  IQ-5!?’  -  273)2  +  1.4454  x  10_7(T  -  273)3. 

(54) 


dr 


-—  (  kG—  -  segp\C\v\T  )  =  hglhVc(RTCv  -  Ps at),  (48) 


dx 


From  these  relationships  we  derive  Cv  as  follows: 

RTcPsat  c  _ 

Cv  =  — 

Pc 


where  Icq  is  the  effective  thermal  conductivity  of  the  GDL. 


(55) 
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where  C  is  the  total  gas  molar  concentration  in  the  cathode  chan¬ 
nel,  calculated  from  the  cathode-channel  pressure: 


(56) 


The  dissolved- water  concentration  at  the  membrane/ ACL  in¬ 
terface  is  assumed  to  be  in  equilibrium  with  the  vapour-phase 
concentration  if  s^Lj)  <  s*,  and  to  be  in  equilibrium  with  water 
if  s(L2)  >  s*.  The  water  activity  at  x  —  Lj,  «w.a,  is  approxi¬ 
mated  by  the  anode-channel  value  (since  water  is  not  produced 
in  the  ACL),  yielding  the  following  expression  for  Cd- 

A.  =  Xd  —  0.3  T  10.8flWja  —  16u^a  -f-  14.1u^,a, 

Cd(i3)  = - — - .  (57) 

1  +  0.0126Aa 


In  the  cathode  channel  the  temperature  is  prescribed.  The 
flux  of  heat  at  x  =  i3  is  approximated  using  the  anode  channel 
temperature  Td.  We  assume  that  the  heat  flux  is  constant  across 
the  anode  GDL  and  ACL,  yielding: 


T(0)  =  Tc, 


dT 

dx 


( LZ)  = 


kc 

Lg 


(Ta  -  T(L3))  . 


(58) 


We  are  assuming  the  that  the  anode  gas  diffusion  layer  is 
identical  to  that  on  the  cathode  side.  This  is  not  necessarily  so 
and  this  assumption  can  easily  be  relaxed. 

The  final  boundary  condition  is  that  for  the  liquid  water  at 
the  interface  between  the  cathode  channel  and  the  GDL.  It  is 
common  in  modelling  studies,  for  example  [3,5-7],  to  assume 
either  a  zero  saturation  or  zero  liquid-water  flux  at  this  location, 
or  along  portions  of  the  channel/GDL  interface  in  two  dimen¬ 
sions.  In  [47],  Weng  and  Wang  specify  the  saturation  at  the 
GDL/channel  interface,  which,  as  they  state,  is  likely  to  depend 
sensitively  on  the  gas  flow  in  the  channel,  current  density  and 
wettability.  According  to  their  results,  high  levels  of  saturation 
are  possible  only  if  the  prescribed  value  at  the  boundary  is  high 
or  the  GDL  permeability  is  unrealistically  small.  In  our  model, 
the  counterpart  to  this  boundary  condition  depends  on  the  gas 
flow  and  wettability,  and  is  written  as  a  flux.  We  show  that  high 
saturation  levels  are  possible  within  a  broad  range  of  realistic 
parameter  space. 

The  accumulation  and  removal  of  liquid  water  along  the  chan¬ 
nel  is  a  complicated  process.  Water  is  expelled  from  the  GDL 
through  preferential  openings  and  forms  droplets  attached  to  the 
surface.  These  droplets  can  grow  or  coalesce  to  a  form  larger 
droplets  comparable  in  size  to  the  channel  dimensions,  which 
results  is  the  formation  of  a  liquid  film,  [48].  Ultimately,  the  film 
is  wicked  along  the  channel  walls  toward  the  exit.  The  greater 
the  flow  velocity  in  the  channel,  the  more  effective  the  removal 
of  liquid  water.  A  detailed  model  of  this  process,  accounting  for 
the  surface  properties  of  the  GDL,  phase  change,  surface  ten¬ 
sion,  two  dimensional  gas  flow  and  hydrophilicity  of  the  chan¬ 
nel  walls,  is  clearly  a  challenge  in  its  own  right  that  will  not  be 
attempted  here.  Instead,  we  approximate  the  process  with  the 
following  steady-state  flux  condition  at  x  =  0: 


(0)  -  G.y(0)  =  0, 

dx 


(59) 


where  =  0  corresponds  to  zero  water  removal.  This  form  is 
motivated  in  two  ways.  Firstly,  there  should  be  no  removal  at  zero 
saturation.  Secondly,  the  parameter  £2  can  be  seen  naturally  as 
(in  an  average  sense)  the  reciprocal  of  a  water  film  thickness  on 
the  surface  of  the  GDL.  When  this  thickness  approaches  zero,  the 
saturation  is  physically  zero,  which  is  expressed  by  1/  £2  — >  0 
in  Eq.  (59).  The  thickness,  and  therefore  £2,  is  likely  to  depend 
sensitively  on  the  flow  rate  in  the  channel,  [48]  (also  by  analogy 
with  the  theory  of  Ekman  transport  in  ocean  modelling). 

4.  Experimental  details 

The  fuel  cell  was  prepared  with  carbon  fibre  paper  (GDL), 
carbon  supported  Platinum  catalyst  and  Nafion  membrane,  with 
a  50  cm2  anode  geometric  area  and  a  masked  cathode  with  a 
4  cm2  geometric  area.  The  MEA  and  cell  hardware  were  de¬ 
signed  to  minimize  any  temperature,  pressure  and  concentration 
gradients  in  the  cell.  The  small  cathode  geometric  area  ensures 
minimal  heat  generation  and  concentration  change  along  the 
length  of  the  flow-field  channels.  To  further  ensure  that  there 
was  no  concentration  change  along  the  length  of  the  flow  field, 
we  used  a  stoichiometry  of  greater  than  60  (compared  to  a  nor¬ 
mal  operating  value  of  approximately  2).  Thus,  the  design  was 
such  that  variations  were  predominantly  through  the  MEA,  and 
therefore  one-dimensional.  Note  also  that  the  flow  rate  was  the 
same  in  all  cases. 

Using  a  load  bank,  a  series  of  polarization  curves  were  pro¬ 
duced  under  varying  conditions  to  demonstrate  the  change  in 
the  polarization  losses  relative  to  oxygen  concentration,  water 
activity  and  temperature. 

The  oxygen  concentration  was  varied  between  3%  and  5%, 
balanced  in  nitrogen.  A  separate  investigation  was  performed 
at  the  same  concentrations  with  oxygen  balanced  with  helium. 
These  results  are  not  discussed  here. 

Three  different  water  activity  conditions  were  tested  at  a  cell 
temperature  of  70  °C:  aw  =  0.6  anode /aw  =  0.7  cathode,  aw  = 
0.7  anode/«w  =  0.8  cathode,  and  aw  =  1  on  both  anode  and 
cathode.  Three  different  cell  temperatures  were  tested  at  aw  =  1 
to  validate  the  effect  of  temperature  on  performance:  70, 45  and 
25  °C.  Only  results  pertaining  to  the  first  are  discussed  here. 

5.  Numerical  results 

5.1.  Finite  element  solution  strategy 

We  first  outline  the  solution  strategy  and  the  manner  in  which 
the  results  were  matched  to  the  experimental  data  (presented  in 
the  next  section).  The  governing  equations  and  boundary  condi¬ 
tions  laid  out  above  were  solved  using  the  finite-element  method 
implemented  in  FEMLAB.  The  discretization  of  the  equations 
was  achieved  on  a  uniform  grid  using  quartic  Lagrange  poly¬ 
nomials,  allowing  the  number  of  grid  points  to  be  kept  small 
(typically  32  or  64).  The  switch  functions  were  substituted  with 
hyperbolic  tanh  functions  to  smooth  the  discontinuities,  a  stan¬ 
dard  procedure.  Adaptivity,  a  feature  of  FEMLAB,  was  not  re¬ 
quired  in  the  calculations,  though  was  occasionally  used  to  ob¬ 
tain  smoother  plots. 
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Fig.  4.  Electrolyte  potential  and  oxygen  concentration  profiles  corresponding  to  aw  =  1  in  Fig.  7,  as  the  cell  potential  is  decreased. 


In  order  to  compare  the  modelling  with  the  experimental  re¬ 
sults,  all  parameter  values  were  fixed  as  specified  below  with  the 
exception  of  the  exchange  current  density,  io2,c-  The  latter  was 
then  adjusted  to  match  the  experimental  data  displayed  in  Figs.  6 
and  7.  This  fitting  was  required  only  once,  i.e.,  z'o2,c  =  40  A  m-2 
in  both  figures.  The  motivation  for  this  fitting  procedure  stems 
from  the  lack  of  a  consistent  value  for  i o,c  in  the  literature 
(426Am-2  in  [41]  and  1.5  x  10-4  Am“2  in  [26]).  Similar  dis¬ 
crepancies  exist  in  the  values  of  the  transfer  coefficients  and 
co2,ref-  The  difficulties  associated  with  measuring  or  estimat¬ 
ing  these  quantities  are  numerous.  For  example,  z'o2,c  changes 
from  one  cathode  structure  to  another  (active  surface  area  of 
platinum),  as  well  as  with  temperature.  The  value  we  derive  for 
z'o2,c  falls  well  within  the  broad  range  of  values  that  appear  in  the 
literature.  A  more  detailed  discussion  of  the  issues  surrounding 
the  measurement  of  these  quantities  can  be  found  in  [54]  and 
[55], 

It  should  be  mentioned  that  Q  was  fixed  at  0.075  m~ 1 .  Since  a 
value  for  £2  is  difficult  to  determine  theoretically,  we  later  present 
a  parametric  study  of  its  effect.  The  flow  rate  in  the  cathode 


channel  is  assumed  to  be  constant  (as  in  the  experiments),  which 
justifies  the  use  of  a  constant  in  all  calculations. 

5.2.  Experimental  validation 

For  the  modelling  results  of  this  section,  the  parameters 
values  that  were  used  are  shown  in  Table  1,  unless  otherwise 
stated. 

We  begin  with  the  set  of  results  demonstrated  in  Figs.  4  and 
5,  the  first  of  which  shows  profiles  of  electrolyte  potential  and 
oxygen  concentration  as  the  cell  voltage  (approximated  by  the 
carbon  potential  Uc)  is  decreased  from  0.95  to  0. 1  V.  In  the  cal¬ 
culations,  both  channels  are  at  aw  =  1  and  70  °C,  and  the  oxygen 
concentration  in  the  cathode  channel  is  5%  of  the  total.  The  sec¬ 
ond  figure  shows  the  saturation  levels  for  both  the  case  above 
and  with  «w  =  0.7  anode/aw  =  0.8  cathode,  indicating  the  ef¬ 
fect  that  the  water  activity  in  the  channels  has  on  the  flooding  of 
the  cathode:  that  saturation  levels  can  be  high  at  low  cell  voltage 
and  with  high  water  activity  in  the  channels.  We  return  to  this 
point  later. 


0  50  100  150  200 

distance  pm 


Fig.  5.  Saturation  levels  for  decreasing  cell  voltage  corresponding  to  «w  =  1  (left)  and  aw  =  0.7  anode/aw  =  0.8  cathode,  the  polarization  curves  for  which  are 
given  in  Fig.  7. 
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Fig.  6.  Comparison  of  the  modelling  and  experimental  polarization  curves  at  5 
and  3%  oxygen  concentration  (see  Table  1  for  a  full  list  of  parameter  values). 


To  demonstrate  the  comparison  between  the  modelling  re¬ 
sults  and  those  of  the  experimental  procedure  described  in  the 
previous  section,  we  refer  to  Figs.  6  and  7.  The  first  of  these  fig¬ 
ures  shows  the  polarization  curves  at  two  different  oxygen  con¬ 
centrations,  with  aw  =  0.7  anode /«w  =  0.8  cathode.  The  sec¬ 
ond  figure  shows  polarization  curves  at  two  different  values  of 
water  activity  in  the  channels,  with  a  5%  oxygen  concentration. 
In  both  of  these  figures,  the  conditions  not  specified  are  as  in 
Table  1,  except  that  the  thickness  of  the  CCL  is  25  p,m. 

The  agreement  between  the  modelling  results  and  the  exper¬ 
imental  data  is  excellent.  The  second  figure  clearly  underlines 
the  ability  of  the  model  to  capture  the  impact  of  water  on  both 
the  membrane  and  CCL/GDL.  The  modelling  and  experimental 
results  strongly  suggest  that  the  increased  hydration  of  the  mem¬ 
brane  at  high  water  activity  leads  to  improved  performance  at 
high  and  mid-range  cell  voltages.  However,  as  the  current  den¬ 


Fig.  7.  Comparison  of  the  modelling  and  experimental  polarization  curves  for 
two  aw  configurations  at  5%  concentration  of  oxygen  (see  Table  1  for  a  full  list 
of  parameter  values). 


sity  increases  so  does  the  water  generated  in  the  CCL,  which 
leads  to  significant  oxygen  diffusion  limitation  from  flooding, 
and  ultimately  to  a  better  performance  at  low  water  activity. 
These  features  are  very  well  captured  by  the  model.  Below  we 
say  more  on  the  effects  of  changes  in  the  channel  water  activity, 
particularly  when  coupled  to  changes  in  the  channel  tempera¬ 
tures. 


5.3.  Effects  of  GDL  and  CCL  structures 


The  magnitudes  of  the  permeabilities,  kc  and  ks,  and  capillary 
diffusion  coefficients,  erc  and  erg,  are  subject  to  a  great  deal  of 
uncertainty.  The  effect  that  changes  in  these  quantities  have  on 
the  flooding  behaviour,  and  therefore  performance,  are  not  well 
understood.  In  this  subsection  make  an  attempt  to  quantify  the 
liquid-water  effects  as  these  quantities  are  varied.  The  benefit 
of  such  information  is  clearly  in  predicting  the  performance  of 
different  media,  namely  in  the  GDL  and  CCL,  whose  proper¬ 
ties  can  be  classified  in  relation  to  contact  angle  (hydrophobic- 
ity/hydrophilicity),  permeability  and  porosity.  Let  us  return  to 
the  definitions  (21),  (22)  and  (50)  and  now  define  the  quantities: 


KC(TC  .  KgOg 

/X  [1 


(60) 


Changes  in  these  quantities  can  arise  from  changes  in  several 
underlying  parameters,  namely  the  absolute  permeabilities,  kc 
and  Kg,  the  porosities,  ep  and  eg,  and  the  contact  angles  6C  and 
0g.  The  permeabilities  can  be  decomposed  using  the  Kozeny- 
Carman  relation,  [32], 
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where  <r/p  c  and  dPig  are  the  mean  pore  diameters  in  the  CCL 
and  GDL,  respectively.  Kc  and  Kg  are  the  so-called  Kozeny 
constants.  Ultimately 
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so  that  dg,  increases  with  increasing  porosity,  contact  angle  (in 
the  hydrophobic  range)  and  mean  pore  diameter  of  the  GDL.  The 
same  is  true  for  dc  except  that  a  decrease  in  the  contact  angle, 
Of  leads  to  an  increase  in  dc.  Thus,  the  results  obtained  from 
varying  dg  and  dc  contain  a  quite  general  degree  of  information. 
For  typical  values  (those  used  by  Mazumder,  [4]),  dg  is  of  the 
order  IIP4- Iff  .  In  the  CCL,  where  the  pores  are  much  smaller 
( 0(  1  0_7m)),  dc  is  likely  to  be  one  or  two  orders  of  magnitude 
smaller  than  dg. 

Fig.  8  demonstrates  the  effect  on  performance  and  liquid- 
water  accumulation  of  changes  in  dg  and  dc.  In  these  calcula¬ 
tions,  both  channels  are  at  aw  =  1  and  other  parameters,  unless 
otherwise  specified,  are  given  in  Table  1 .  It  can  be  seen  from  the 
first  of  these  figures  that  increases  in  dg  and  dc  lead  to  decreases 
in  the  maximum  saturation  and  increases  in  the  minimum  satura¬ 
tion  (minimum  at  jr  =  0  and  maximum  at  x  —  Lf),  accompanied 
by  increases  in  the  current  density,  measured  at  0. 1  V.  Changes 
in  dg,  at  fixed  dc,  appear  to  have  a  more  dramatic  effect  than 
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Fig.  8.  The  effect  on  the  maximum  saturation  and  current  density  at  0.1  and  0.6  V  of  changes  in  the  contact  angles,  absolute  permeabilities  and  porosities  of  the  CCL 
and  GDL.  dg  and  dc  are  defined  in  (60)  and  both  channels  are  at  aw  =  1 . 


changes  in  dc,  at  fixed  dg,  particularly  on  the  minimum  satura¬ 
tion.  It  can  further  be  seen  that  for  a  fixed  dc  there  is  a  lower 
limit  to  the  saturation  and  an  upper  limit  to  the  current  density 
as  dg  is  increased.  The  difference  in  the  lower  limit  of  maximum 
saturation  as  dc  is  increased  is  quite  appreciable,  but  the  effect 
on  the  current  density,  measured  at  both  0. 1  and  0.6  V,  is  neg¬ 
ligible.  It  appears  then  that  changes  in  the  contact  angle,  pore 
radius  and  porosity  in  the  GDL  have  a  far  greater  impact  than 
changes  in  the  same  quantities  in  the  CCL  (possibly  because  of 
the  far  greater  thickness  of  the  GDL).  In  fact  they  can  alter  the 
accumulation  of  liquid  water  and  consequently  the  performance 
quite  dramatically. 

We  note  finally  that  in  the  range  —8.5  <  logi0  dg  <  —6,  the 
current  density  measured  at  0.6  V  does  not  decrease  as  rapidly 
as  the  current  density  at  0. 1  V.  This  appears  once  more  to  be  the 
consequence  of  the  competition  between  flooding  and  hydration 
of  the  membrane. 


5.4.  Effects  of  GDL  and  membrane  thickness 

The  effects  of  membrane  and  GDL  thickness  changes  can 
be  inferred  from  Fig.  9.  For  these  calculations  both  channels 
are  at  aw  —  1  and  other  parameters  are  given  in  Table  1,  unless 
otherwise  stated.  The  left-hand  figure  demonstrates  the  polar¬ 
ization  curves  for  an  increasing  range  of  membrane  thickness, 
with  a  GDL  thickness  of  200  pun.  It  can  be  seen  from  the  right- 
hand  figure  that  decreasing  Lm  improves  the  performance  by 
increasing  the  slope  of  the  polarization  curve.  Quite  substantial 
improvements  are  possible.  The  right-hand  figure  demonstrates 
that  decreasing  the  GDL  thickness  improves  performance.  Inter¬ 
estingly,  this  is  despite  a  simultaneous  increase  in  the  maximum 
saturation.  This  implies  that  improving  the  oxygen  supply  to  the 
CCL,  by  reducing  the  mass-transport  resistance  from  the  GDL, 
can  to  some  extent  mitigate  the  effects  of  flooding,  another  ex¬ 
ample  of  the  competing  effects  in  the  MEA. 
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Fig.  9.  The  left-hand  figure  shows  the  polarization  curves  for  a  decreasing  range  of  membrane  thickness,  with  a  GDL  thickness  of  200  fjim  and  both  channels  at 
aw  =  1 .  The  right-hand  figure  shows  the  effect  of  changes  in  the  GDL  thickness  on  both  the  maximum  saturation  and  current  density  measured  at  a  cell  voltage  of 
0.1  V.  In  this  plot  Lm  is  given  in  Table  1. 
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Fig.  10.  The  effect  on  the  performance  of  changes  in  the  channel  water  activity,  when  maintained  at  the  same  value  in  both  channels  (left  figure)  and  at  different 
values  (right  figure). 


5.5.  Effects  of  water  activity 

In  Fig.  10  we  plot  the  polarization  curves  for  several  values 
of  water  activity  in  the  channels.  The  left-hand  figure  demon¬ 
strates  that  changes  in  water  activity,  with  the  same  values  in 
both  channels,  have  a  convoluted  effect  on  the  performance.  At 
fully  humidified  conditions  the  effects  of  flooding  are  strongly 
felt  at  low  cell  voltages.  However,  the  performance  at  high 
and  mid-range  voltages  is  good.  At  lower  humidification  lev¬ 
els  down  to  approximately  aw  =  0.8  the  performance  improves 
at  low  cell  voltages,  with  minor  deterioration  at  higher  voltages. 
However,  when  the  water  activity  is  reduced  further,  the  per¬ 
formance  deteriorates,  particularly  at  high  and  mid-range  cell 
voltages  but  also  now  at  lower  cell  voltages  (compared  to  aw  = 
0.8). 

Maintaining  the  cathode-channel  water  activity  and  varying 
that  for  the  anode  channel  leads  to  the  polarization  curves  in 
the  right-hand  plot  of  Fig.  10.  The  best  performance  is  seen  at 
aw, a  =  1,  with  only  a  minor  deterioration  down  to  aw  a  =  0.8. 
Once  below  this  value,  the  increased  resistivity  of  the  mem¬ 
brane  degrades  the  performance  substantially,  particularly  at 
mid-range  cell  voltages. 

5.6.  Combined  effects  of  temperature  and  water  activity 

Fig.  1 1  demonstrates  the  effects  of  simultaneous  changes  in 
the  operating  temperatures  and  water  activities,  at  constant  pres¬ 
sure  in  the  channels.  It  appears  that  at  as  we  increase  temperature, 
changes  in  water  activity  are  less  strongly  felt.  At  T  =  60  °C, 
significant  improvements  in  performance  are  possible  at  mid¬ 
range  and  high  cell  voltages  if  the  water  activity  of  the  channels  is 
high,  as  already  seen  in  the  previous  subsection.  At  T  =  120  °C, 
while  improvements  are  still  possible,  they  are  not  as  dramatic. 
One  implication  of  this  result  is  that  in  some  measure  it  is  possi¬ 
ble  to  mitigate  the  effects  of  flooding  (dramatic  deterioration  in 
performance  at  low  cell  voltage)  by  increasing  the  temperature 
in  the  channels,  and  therefore  the  average  temperature  of  the 


system.  In  Fig.  12  we  see  the  reason  for  this.  At  T  =  60  °C  the 
saturation  levels  are  high  at  low  cell  voltage,  causing  the  drop 
in  performance  seen  in  Fig.  1 1 .  However,  due  to  decreased  con¬ 
densation,  at  T  =  120  °C  even  at  low  cell  voltage  the  saturations 
levels  are  too  low  to  lead  to  flooding. 

We  note  that  in  Fig.  11,  the  gas  pressure  in  the  channel  is 
constant,  so  an  increase  in  temperature  leads  to  a  decrease  in  the 
molar  density  of  the  gas,  and  therefore  in  total  oxygen  avail¬ 
able.  More  importantly,  an  increase  in  temperature  causes  a 
substantial  increase  in  the  saturation  vapour  pressure  at  high 
temperatures.  At  fixed  water  activity,  there  are  correspondingly 
large  increases  in  the  vapour  density.  Fig.  13  shows  the  polar¬ 
ization  curves  at  aw  =  1  as  the  temperature  is  increased  from 
60  to  120  °C,  while  the  gas  pressure,  Pc,  and  the  concentra¬ 
tion  of  oxygen  in  the  cathode  channel,  Cp  are  kept  constant  at 
30psig  (206.842 kPa)  and  3.75 mol m-3,  respectively.  We  see 
that  the  current  density  increases  dramatically  at  low  cell  volt¬ 
age  because  of  the  reduction  in  the  saturation  levels.  At  higher 
cell  voltages  the  performance  deteriorates,  but  not  substantially, 
possibly  because  of  the  increase  in  the  vapour  density  and  its 
corresponding  effect  on  the  membrane  conductivity. 

5.7.  Effects  of  water  removal  from  channel 

To  see  the  effect  of  the  removal  of  liquid  water  from  the 
cathode  channel  we  refer  to  Fig.  14.  It  is  evident  that  by  in¬ 
creasing  the  removal  rate  (equivalent  to  increasing  £2)  perfor¬ 
mance  improves  at  0. 1  V.  The  improvement  in  performance  is 
accompanied  by  a  decrease  in  the  maximum  saturation  (over 
GDL  and  CCL),  and  there  is  an  inverse  relationship  between  the 
flooding  level  and  the  current  density  at  this  voltage.  Moreover, 
the  current  density  measured  at  0.6  V  does  not  decrease  signifi¬ 
cantly  over  the  same  range  of  12.  In  other  words,  increasing  the 
removal  of  liquid  water  from  the  channel  improves  the  perfor¬ 
mance  at  low  cell  voltages  without  significantly  affecting  the 
performance  at  high  cell  voltages.  This  point  is  underlined  by 
Fig.  15. 
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Fig.  11.  The  effect  on  the  performance  of  combined  changes  in  the  channel  water  activity  and  temperature  (both  channels  at  the  same  water  activity). 


6.  Discussion 

Capturing  the  effects  of  (particularly  liquid)  water  is  a  non¬ 
trivial  issue  in  PEMFC  modelling.  In  this  paper  we  have  provided 
the  outline  of  a  model  that  incorporates  many  of  the  difficult 
aspects  associated  with  these  effects.  In  doing  so,  we  found  it 
necessary  to  explicitly  include  many  of  the  MEA  features  that 
are  often  subject  to  approximation,  such  as  a  detailed  catalyst 
layer  model,  temperature  variations  and  the  multi-phase  nature 
of  the  problem. 

We  have  demonstrated  through  comparison  with  experi¬ 
ment  that  the  model  we  have  developed  qualitatively  captures 
the  correct  behaviour  with  respect  to  flooding  and  membrane 
hydration  (see  Fig.  7).  This  result,  which  appears  to  be  ab¬ 
sent  in  previous  modelling  studies,  underlines  the  competi¬ 
tion  between  hydration  of  the  membrane  and  flooding,  and  is 
clearly  an  effect  which  has  to  be  captured  in  any  model  that 
is  aimed  at  predicting  or  aiding  water  management.  To  show 
the  benefits  of  our  modelling  approach  we  performed  several 
studies  focussed  on  water  management,  particularly  with  re¬ 
gard  to  the  effects  that  micro-structure  and  operating  condi¬ 
tions  have  on  performance.  We  now  highlight  a  few  of  our 
findings: 


(1)  From  the  results  depicted  in  Fig.  8  and  the  relation¬ 
ships  (60)-(62),  it  appears  that  the  most  influential  micro- 
structural  property  to  affect  saturation  levels  is  the  pore  di¬ 
ameter,  d? g  for  the  GDL  and  dp  c  for  the  CCL.  Provided 
that  they  do  not  approach  extreme  values  (#£,  0§  — »■  90°), 
the  contact  angles  play  only  a  minor  role  since  the  maxi¬ 
mum  saturation  is  quite  insensitive  to  changes  in  dc  and  dg 
that  are  not  at  least  an  order  of  magnitude.  Nevertheless, 
combined  changes  in  the  microstructural  properties  of  the 
GDL  can  be  employed  to  improve  performance  quite  sub¬ 
stantially. 

(2)  Reducing  the  thickness  of  the  membrane  improves  perfor¬ 
mance.  However,  this  practice  must  be  treated  with  cau¬ 
tion.  If  the  thickness  is  reduced  by  too  much,  oxygen  and 
hydrogen  cross-over  will  also  increase  and  will  lead  to  se¬ 
vere  degradation.  Likewise,  decreasing  the  GDL  thickness 
improves  performance,  with  the  increased  oxygen  supply 
outweighing  the  increased  saturation  levels. 

(3)  From  the  results  of  Section  5.5,  it  is  clear  that  without  either 
changes  in  the  other  conditions  or  external  humidification 
strategies,  the  resistance  of  the  membrane  is  prohibitively 
high  if  the  water  activities  in  the  channels  are  below  a  critical 
set  of  values  (around  0.7-0. 8).  At  water  activities  approach- 
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Fig.  12.  Saturation  profiles  corresponding  to  the  calculations  at  aw  =  1  in  Fig.  11. 


ing  unity  in  both  channels,  the  flooding  becomes  an  issue  if 
the  cell  is  operated  at  low  voltage.  However,  at  mid-range 
voltages,  above  approximately  0.4  V,  the  performance  is  op¬ 
timized  at  ctw  =  1 . 

(4)  If  the  cell  is  operated  at  water  activities  in  the  channel  that 
are  close  to  unity,  temperature  can  be  used  in  some  measure 
to  offset  the  effects  of  flooding  by  lowering  the  condensa¬ 
tion  rate  and  therefore  the  saturation  levels  (see  Figs.  1 1 
and  12).  From  Fig.  13  we  see  that  performance  at  low  cell 
voltages  can  therefore  be  enhanced  significantly,  while  the 
corresponding  deterioration  in  performance  at  higher  cell 
voltages  is  minor.  The  latter  result  is  possibly  due  to  the  in¬ 
creased  vapour  molar  concentration  at  higher  temperature 
and  the  corresponding  improvement  in  membrane  conduc¬ 
tivity. 

(5)  Fig.  14  suggests  that  at  high  water  activity  significant  per¬ 
formance  improvements  can  be  made  by  increasing  the  re¬ 
moval  rate  of  liquid  water  from  the  cathode  channel,  which 
although  not  solely  controlled  by,  is  directly  linked  to  the 
flow  rate  of  gas  down  the  channel.  Under  these  conditions, 
if  no  liquid  water  is  removed  (!T  =  0)  flooding  severely  im¬ 
pacts  the  current  density  at  low  cell  voltage.  The  greater  the 
removal  of  the  liquid  water  the  higher  the  current  density  at 


Fig.  13.  The  effect  of  temperature  increases  in  the  channels  (both  at  the  same 
temperature)  with  aw  =  1  in  both  channels.  In  these  calculations,  the  oxygen 
concentration  was  kept  constant  as  temperature  was  varied. 
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Fig.  14.  The  effect  of  the  removal  constant  £2,  defined  in  (59),  on  the  maximum 
saturation  and  the  current  density  measured  at  0.1  and  0.6  V.  Both  channels  are 
at  aw  =  1 . 


low  cell  voltage,  which  eventually  plateaus  as  Q  increases. 
Furthermore,  the  performance  at  higher  cell  voltage  is  not 
dramatically  affected  as  the  removal  increases  (see  Fig.  15). 


Though  we  have  concentrated  on  water  effects  in  this  paper 
there  are  several  other  aspects  to  the  model  that  can  be  exploited. 
For  example,  a  parametric  study  with  respect  to  the  microstruc- 
tural  properties  of  the  CCL  and  conditions  under  which  proton 
migration  is  limiting.  Our  purpose  has  been  to  demonstrate  the 
ability  of  the  model  to  capture  complex  phenomena.  There  are 
several  extensions  to  the  model  that  we  are  currently  pursuing. 
In  order  to  study  effects  such  as  hydrogen  and  oxygen  (as  well 
as  water)  crossover,  the  anode  catalyst  and  gas-diffusion  layers 
must  be  included.  In  addition,  extension  to  two  dimensions  will 
provide  details  of  the  non-uniformities  down  the  channels. 


Fig.  15.  The  effect  on  the  polarization  curve  of  changes  in  £2,  the  removal  rate 
constant. 


Acknowledgements 

The  authors  are  grateful  to  MITACS  and  Ballard  Power  Sys¬ 
tems  for  their  financial  support  and  experimental  facilities.  We 
would  particularly  like  to  acknowledge  Dr.  Jean  St.  Pierre  (for¬ 
merly  of  Ballard  Power  Systems),  Professor  Brian  Wetton  and 
Dr.  John  Stockie  (MITACS)  for  their  involvement  in  the  many 
fruitful  discussions  leading  to  this  paper. 


References 

[1]  J.  Cao,  N.  Djilali,  Computational  simulation  of  water  transport  in  PEM 
fuel  cells  using  an  improved  membrane  model.  Proceedings  of  the  10th 
Canadian  Hydrogen  Conference  (2000)  447^-56. 

[2]  W.  He,  J.  Yi,  T.V.  Nguyen,  Two-phase  flow  model  of  the  cathode  of  PEM 
fuel  cells  using  interdigitated  flow  fields,  AIChE  J.  46  (10)  (2000)  2053- 
2064. 

[3]  G.  Lin,  W.  He,  T.  Van  Nguyen,  Modelling  liquid  water  effects  in  the  gas  dif¬ 
fusion  and  catalyst  layers  of  the  cathode  of  a  PEM  fuel  cell,  J.  Electrochem. 
Soc.  151  (11)  (2004)  A1999-A2006. 

[4]  S.  Mazumder,  J.V.  Cole,  Rigorous  3D  mathematical  modelling  of  PEM 
fuel  cells:  2  model  predictions  with  liquid  water  transport,  J.  Electrochem. 
Soc.  150  (11)  (2004)  A1510-A1517. 

[5]  D.  Natarajan,  T.V.  Nguyen,  A  two-dimensional,  two-phase,  multicompo¬ 
nent,  transient  model  for  the  cathode  of  a  proton  exchange  membrane  fuel 
cell  using  conventional  gas  distributors,  J.  Electrochem.  Soc.  148  (12) 
(2001)  A1324-A1335. 

[6]  N.P.  Siegel,  M.W.  Ellis,  D.J.  Nelson,  M.R.  von  Spakovsky,  Single  domain 
PEMFC  model  based  on  agglomerate  catalyst  geometry,  J.  Power  Sources 
115  (1)  (2003)  81-89. 

[7]  N.P.  Siegel,  M.W.  Ellis,  D.J.  Nelson,  M.R.  von  Spakovsky,  A  two- 
dimensional  computational  model  of  a  PEMFC  with  liquid  water  transport, 
J.  Power  Sources  128  (1)  (2004)  173-184. 

[8]  L.B.  Wang,  N.I.  Wakayama,  T.  Okada,  Numerical  simulation  of  enhance¬ 
ment  of  mass  transfer  in  the  cathode  electrode  of  a  PEM  fuel  cell  by  magnet 
particles  deposited  in  the  cathode-side  catalyst  layer,  Chem.  Eng.  Sci.  60 
(2005)  4453-4467. 

[9]  Z.H.  Wang,  C-Y.  Wang,  K.S.  Chen,  Two-phase  flow  and  transport  in  the 
air  cathode  of  proton  exchange  membrane  fuel  cells,  J.  Power  Sources  94 
(2001)  40-50. 

[10]  S.  Lister,  G.  McLean,  Review:  PEM  fuel  cell  electrodes,  J.  Power  Sources 
130  (2004)  61-76. 

[11]  F.  Liu,  B.  Yi,  D.  Xing,  J.  Yu,  Z.  Hou,  Y.  Fu,  Development  of  novel  self- 
humidifying  composite  membranes  for  fuel  cells,  J.  Power  Sources  124 
(2003)  81-89. 

[12]  M.  Uchida,  Y.  Aoyama,  N.  Eda,  O.  Akira,  Investigation  of  the  microstruc¬ 
ture  in  the  catalyst  layer  and  effects  of  both  perfluorosulfonate  ionomer 
and  PTFE-loaded  carbon  on  the  catalyst  layer  of  polymer  electrolyte  fuel 
cells,  J.  Electrochem.  Soc.  142  (12)  (1995)  A4143-A4149. 

[13]  M.  Uchida,  Y.  Fukuota,  Y.  Sugawara,  H.  Ohara,  O.  Akira,  Improved 
preparation  process  of  very-low-platinum-loading  electrodes  for  poly¬ 
mer  electrolyte  fuel  cells,  J.  Electrochem.  Soc.  145  (11)  (1998)  A3708- 
A3713. 

[14]  A.Z.  Weber,  J.  Newman,  Modelling  transport  in  polymer-electrolyte  fuel 
cells,  Chem.  Rev.  104  (2004)  4679^1726. 

[15]  K.  Broka,  P.  Ekdunge,  Modelling  the  PEM  fuel  cell  cathode,  J.  Appl.  Elec¬ 
trochem.  27  (1997) 281-289. 

[16]  Q.  Guo,  R.  White,  A  steady-state  impedance  model  for  a  PEMFC  cathode, 
J.  Electrochem.  Soc.  151  (4)  (2004)  E133-E149. 

[17]  F.  Jaouen,  G.  Lindbergh,  G.  Sundholm,  Investigation  of  mass-transport 
limitations  in  the  solid  polymer  fuel  cell  cathode:  I.  mathematical  model, 
J.  Electrochem.  Soc.  149  (4)  (2002)  A437-A447. 

[18]  S.  Mazumder,  J.V.  Cole,  Rigorous  3D  mathematical  modelling  of  PEM  fuel 
cells:  1  model  predictions  without  liquid  water  transport,  J.  Electrochem. 
Soc.  150  (11)  (2004)  A1503-A1509. 


1268 


A.A.  Shah  et  al.  /  Journal  of  Power  Sources  160  (2006)  1251-1268 


[19]  W.  Sun,  B.A.  Peppley,  K.  Karan,  An  improved  two-dimensional  agglom¬ 
erate  cathode  model  to  study  the  influence  of  catalyst  layer  structural  pa¬ 
rameters,  Electrochim.  Acta  50  (2005)  3359-3374. 

[20]  K.  Yin,  Parametric  study  of  proton-exchange  membrane  fuel-cell  cathode 
using  an  agglomerate  model,  J.  Electrochem.  Soc.  152  (3)  (2005)  A583- 
A593. 

[21]  R.  Aris,  The  Mathematical  Theory  of  Diffusion  and  Reaction  in  Permeable 
Catalysts,  vol.  1.  The  Theory  of  the  Steady  State,  Oxford  University  Press, 
1975. 

[22]  M.  Eikerling,  A.A.  Kornyshev,  Modelling  the  performance  of  the  cathode 
catalyst  layer  of  polymer  electrolyte  fuel  cells,  J.  Electroanal.  Chem.  453 
(1998)  89-106. 

[23]  M.  Eikerling,  A.A.  Kornyshev,  Electrochemical  impedance  of  the  cathode 
catalyst  layer  in  polymer  electrolyte  fuel  cells,  Electroanal.  Chem.  475 
(1999)  107-123. 

[24]  D.B.  Genevey,  Transient  model  of  heat,  mass  and  charge  transfer  as  well 
as  electrochemistry  in  the  cathode  catalyst  layer  of  a  PEMFC,  Master’s 
Thesis,  Virginia  Polytechnic  Institute,  2004. 

[25]  H.-K.  Hsuen,  Mechanistic  approach  to  performance  equations  for  cath¬ 
odes  in  polymer  electrolyte  fuel  cells,  J.  Power  Sources  123  (2003)  26- 
36. 

[26]  D.  Song,  Q.  Wang,  Z.  Liu,  T.  Navessin,  M.  Eikerling,  S.  Holdcroft,  Nu¬ 
merical  optimization  study  of  the  catalyst  layer  of  a  PEM  fuel  cell  cathode, 
J.  Power  Sources  126  (2004)  104-1 11. 

[27]  C.Y.  Sukkee  Um,  K.S.  Wang,  Chen,  Computational  fluid  dynamics  mod¬ 
eling  of  proton  exchange  membrane  fuel  cells,  J.  Electrochem.  Soc.  147 
(12)  (2000)  A4485-A4493. 

[28]  Q.  Wang,  M.  Eikerling,  D.  Song,  Z.  Liu,  Structure  and  performance  of 
different  types  of  agglomerates  in  cathode  catalyst  layer  of  PEM  fuel  cells, 
J.  Electroanal.  Chem.  573  (2004)  61-79. 

[29]  H.  Ju,  H.  Meng,  C.-Y.  Wang,  A  single-phase,  non-isothermal  model 
for  PEM  fuel  cells,  Int.  J.  Heat  Mass  Transfer  148  (7)  (2005)  1303 
-1315. 

[30]  H.  Weng,  C.-Y.  Wang,  Model  of  two-phase  flow  and  flooding  dynamics  in 
polymer  electrolyte  fuel  cells,  J.  Electrochem.  Soc.  152  (9)  (2005)  A 1733- 
A1741. 

[31]  C.-Y.  Wang,  P.  Cheng,  Multiphase  flow  and  heat  transfer  in  porous  media, 
J.  Power  Sources  30  (1997)  93-196. 

[32]  J-H.  Nam,  M.  Kaviany,  Effective  diffusivity  and  water- saturation  distribu¬ 
tion  in  single-  and  two-layer  PEMFC  diffusion  medium,  Int.  J.  Heat  Mass 
Trans.  46  (2003)  4595^4611. 

[33]  U.  Pasaogullari,  C-Y.  Wang,  K.S.  Chen,  Two-phase  transport  in  polymer 
electrolyte  fuel  cells  with  bilayer  cathode  diffusion  media,  J.  Electrochem. 
Soc.  152  (8)  (2005)  A1574-A1582. 

[34]  E.  Middelman,  Improved  PEM  fuel  cell  electrodes  by  controlled  self- 
assembly,  Fuel  Cells  Bull.  (2002)  9-12. 

[35]  H.  Mizuhata,  N.  Shin-ichi,  T.  Yamaguchi,  Morphological  control  of 
PEMFC  electrode  by  graft  polymerization  of  polymer  electrolyte  onto 
platinum-supported  carbon  black,  J.  Power  Sources  138  (2004)  25-30. 

[36]  M.  Mathias,  J.  Roth,  J.  Fleming,  W.  Lehnert,  Diffusion  media  materials  and 
characterization,  in:  W.  Vielstich,  A.  Lamm,  H.  Gasteiger  (Eds.),  Handbook 


of  Fuel  Cells — Fundamentals,  Technology  and  Applications,  vol.  3,  John 
Wiley  &  Sons,  2003  (Chapter  46). 

[37]  B.  Bird,  W.  Stewart,  E.  Lightfoot,  Transport  Phenomena,  John  Wiley  and 
Sons,  2002. 

[38]  Z.  Ogumi,  Z.  Takehara,  S.  Yoshizawa,  Gas  permeation  in  spe  method  I. 
oxygen  permeation  through  nation  and  neosepta,  J.  Electrochem.  Soc.  131 
(4)  (1984)  A769-A772. 

[39]  T.E.  Springer,  T.A.  Zawodinski,  S.  Gottesfeld,  Polymer  electrolyte  fuel  cell 
model,  J.  Elecrochem.  Soc.  138  (4)  (1991)  A2334-A2341. 

[40]  T.  Thapman,  S.  Malhotra,  H.  Tang,  R.  Datta,  Modeling  of  conductive  trans¬ 
port  in  proton-exchange  membranes  for  fuel  cells,  J.  Elecrochem.  Soc.  147 
(9)  (2000)  A3242-A3250. 

[41]  P.  Berg,  K.  Promislow,  J.  St.  Pierre,  J.  Stumper,  B.  Wetton,  Water  manage¬ 
ment  in  PEM  fuel  cells,  J.  Electrochem.  Soc.  151  (3)  (2004)  A341-A353. 

[42]  S.  Motupally,  A.J.  Becker,  J.W.  Weidner,  Diffusion  of  water  in  nafion  115 
membranes,  J.  Electrochem.  Soc.  147  (9)  (2000)  A3171-A3177. 

[43]  J.T.  Hinatsu,  M.  Mizuhata,  H.  Takenaka,  Water  uptake  of  perfluorosulfonic 
acid  membranes  from  liquid  water  and  water  vapor,  J.  Elecrochem.  Soc. 
141  (1994)  A1493-A1497. 

[44]  S.  Ge,  X.  Li,  B.  Yi,  I.M.  Hsing,  Absorption,  desorption  and  transport  of 
water  in  polymer  electrolyte  membranes  for  fuel  cells,  J.  Electrochem.  Soc. 
152  (6)  (2005)  A1149-A1157. 

[45]  A.Z.  Weber,  J.  Newman,  Transport  in  polymer-electrolyte  membranes  II. 
mathematical  model,  J.  Elecrochem.  Soc.  151  (2)  (2004)  A311-A325. 

[46]  S.  Mazumder,  A  generalized  phenomenological  model  and  database  for 
the  transport  of  water  and  current  in  polymer  electrolyte  membranes,  J. 
Electrochem.  Soc.  152  (8)  (2005)  A1633-A1644. 

[47]  H.  Weng,  C-Y.  Wang,  Electron  transport  in  pefcs,  J.  Electrochem.  Soc.  151 
(3)  (2004)  A3 5 8- A3 67. 

[48]  X.G.  Yang,  F.Y.  Zhang,  A.L.  Lubawy,  C.Y.  Wang,  Visualization  of  liquid 
water  transport  in  a  PEFC,  Electrochem.  Solid-State  Lett.  7  (11)  (2004) 
A408-A411. 

[49]  M.V.  Williams,  H.R.  Kunz,  J.M.  Fenton,  Influence  of  convection  through 
gas-diffusion  layers  on  limiting  current  in  PEM  fcs  using  a  serpentine  flow 
field,  J.  Elecrochem.  Soc.  151  (10)  (2004)  A1617-A1627. 

[50]  T.  Van  Nguyen,  W.  He,  Interdigitated  flow  field  design,  in:  W.  Vielstich, 
A.  Lamm,  H.  Gasteiger  (Eds.),  Handbook  of  Fuel  Cells — Fundamentals, 
Technology  and  Applications,  vol.  3,  John  Wiley  &  Sons,  2003  (Chapter 
46). 

[51]  R.H.  Perry,  D.W.  Green,  J.O.  Maloney,  Perry’s  Chemical  Engineers  Hand¬ 
book,  McGraw-Hill,  1984. 

[52]  Suk  Won  Cha,  Scaling  effects  of  flow  channels  in  fuel  cells,  PhD  Thesis, 
Stanford  University,  2003. 

[53]  M.J.  Lampinen,  M.  Fomino,  Analysis  of  free  energy  and  entropy  changes 
for  half-cell  reactions,  J.  Electrochem.  Soc.  140  (12)  (1993)  A3537 
-A3546. 

[54]  Q.  Guo,  V.A.  Sethuraman,  R.  White,  Parameter  estimates  for  a  PEMFC 
cathode,  J.  Electrochem.  Soc.  151  (7)  (2004)  A983-A993. 

[55]  B.  Carnes,  N.  Djilal,  Systematic  parameter  estimation  for  PEM  fuel  cell 
models,  J.  Power  Sources  144  (1)  (2005)  83-93. 


